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ABSTRACT 

The magnetorotational instability (MRI) is a crucial mechanism of angular momentum transport in 
a variety of astrophysical accretion disks. In systems accreting at well below the Eddington rate, such 
as the central black hole in the Milky Way (Sgr A*), the rate of Coulomb collisions between particles 
is very small, making the disk evolve essentially as a collisionless plasma. We present a nonlinear 
study of the collisionless MRI using first-principles particle-in-cell (PIC) plasma simulations. In this 
initial study we focus on local two-dimensional (axisymmetric) simulations, deferring more realistic 
three-dimensional simulations to future work. For simulations with net vertical magnetic flux, the 
MRI continuously amplifies the magnetic field, B, until the Alfven velocity, va, is comparable to 
the speed of light, c (independent of the initial value of va/c). This is consistent with the lack of 
saturation of MRI channel modes in analogous axisymmetric MHD simulations. The amplification of 



the magnetic field by the MRI generates a significant pressure anisotropy in the plasma (pj 



> 



P\\, 



where p± and p\\ are the plasma pressures perpendicular and parallel to the local magnetic field). We 
find that this pressure anisotropy in turn excites mirror modes and that the volume averaged pressure 
anisotropy remains near the threshold for mirror mode excitation. Particle energization is due to both 
reconnection and viscous heating associated with the pressure anisotropy. Reconnection produces 
a distinctive power-law component in the energy distribution function of the particles, indicating 
the likelihood of non-thermal ion and electron acceleration in collisionless accretion disks. This has 
important implications for interpreting the observed emission - from the radio to the gamma-rays - 
of systems such as Sgr A*. 

Subject headings: accretion disks, MRI, kinetic plasma effects 
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1. INTRODUCTION 

Accretion disks are ubiquitous in astrophysics and 
play a fundamental role in areas as diverse as planet 
formation, gamma ray bursts (GRBs), and accretion 
onto super massive black holes in the centers of galaxies. 
The accretion of gas in disks requires outward trans- 
port of angular momentum, typically assumed to be 
provided by the magnetoro tational instability (MRI; 
iBalbus fc Hawle^[T^ . [T9^ . The MRI has been widely 
studied using MHD simulations. However, in many cases 
the MHD approach is not directly applicable. When 
the time scale for electron and ion Coulomb collisions 
is longer than the inflow time in the disk, the plasma 
is macroscopically collisionless and MHD breaks down. 
This is the case in radiatively inefficient accretion flow 
(RIAF) models, applicable when the ac cretion rate is 
less than a few percent of Eddington (jNaravan et al.l 
1998). The low rate of Coulomb collisions implies 
that ions and electrons are thermally decoupled, so 
the plasma should be two- temperature. RIAFs are 
ubiquitous, occurring, for example, in t he low-hard state 
of X-ray binaries (e.g. JEsin et al.| [l997). and around the 
central super massive black hole in the Milky Way (Sgr 



A*) and most nearby galaxies. 

The first efforts to understand the MRI in the collision- 
less limit were done using the kinetic MHD approach 
(Quataert et al. 2002; Sharma et al. 2003; Sharma et al. 
2006; see also the closely related work by Balbus 2004 
and Islam & Balbus 2005). These studies highlighted 
the importance of pressure anisotropics with respect 
to the local magnetic field, and incorporated the 
evolution of the pressure parallel and perpendicular to B 
in a fluid model of the plasma. In particular, an increase 
in the magnetic field due to the MRI causes the plasma 
temperature Tj (and pressure, pj) of a given particle 
species j to increase in the direction perpendicular to 
the local magnetic field due to the conservation of the 
magnetic moment, jij, on scales larger than their Larmor 
radius (/ij = p±j/pjB, where pj is the mass density 
of species j). (Throughout this paper the subscript j 
will stand for the particle species, j = i for ions and 
j = e for electrons). This way, on average, the pressures 
perpendicular (p±j) and parallel (p\\j) to the local 
magnetic field must satisfy p±j > p\\j. We expect 
that this anisotropy is ultimately regulated by kinetic 
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micro-instabilities (e.g., ion cyclotron, mirror, firehose, 
electron whistler, etc), as shown by previous detailed 
calculations of the kinetic stability of pl asmas with 
aniso tropic pressure, PIC simulations (e.g., iGary et aTl 
119971 ). and by solar wind observations (Bale et al. 2009). 
The kinetic MHD simulations of Sharm aet al.l (2006, 
2007) modeled the presence of these instabilities by 
setting an upper limit to \T±j/T\\j — 1|. This limit on 
the temperature anisotropy plays a critical role in the 
evolution of the MRI, making the physics much more 
MHD-like than it otherwise would have been. Indeed, 
the satura tion of the MRI is qualitati vely similar to that 
in MHD (jSharma et al.l I2006L 12007ft . One significant 
difference, however, is that the presence of a pressure 
anisotropy leads to an anisotropic pressure stress that 
may be as important for angular momentum transport 
and plasma heating as the magnetic stress. 

In this paper, we study the collisionless MRI using first- 
principle ID and 2D particle-in-cell (PIC) simulations. 
We defer more realistic - but also more computationally 
expensive - 3D calculations to a future paper. In a 
PIC code, the plasma is represented by a collection of 
macro-particles that carry charge and mass, and are 
moved by integration of the Lorentz equations. The 
electromagnetic fields are evolved by solving Maxwell's 
equations on a grid, where the current is calculated 
by adding the contribution of each macro-particle. 
Given its complete treatment of plasmas, the PIC 
approach has the ability to capture the whole dynamics 
of the particles and fields. In particular, the MRI, the 
resulting evolution of the plasma pressure anisotropy, 
the interaction of particles with small-scale kinetic 
instabilities, and particle heating and acceleration are 
all self-consistently captured. 

The reminder of this paper is organized as follows. The 
basic equations and the simulation setup are explained in 
$2] and respectively. In a thorough dispersion rela- 
tion study is performed based on ID simulations, which 
we compare with previous analytic results. We also study 
the nonlinear magnetic field amplification by the MRI 
in ID. In £[5j we explore the non-linear evolution of the 
MRI-driven turbulence using more realistic 2D simula- 
tions. Special attention is paid to the field saturation, 
particle heating, pressure anisotropy evolution, and the 
contribution of the different stress tensor components to 
angular momentum transport. Finally, in $6]we present 
our conclusions. 



acquire extra terms (SchifT 1939): 

V .E = 4np c + ^^--^-VxB, (1) 

c c 

V-B = 0, (2) 
dB -> 

— = -cV x E and (3) 
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where vq = ujq x f and c is the speed of light. 
Our approach is to neglect all of the terms due to the 
non-inertial reference frame in equations [T] and HJ We 
now justify this approximation. In the non-relativistic 
limit (vq <C c), the last two terms on the right hand side 

of Equation [4] are much smaller than cV x B. Thus, when 
<C c, it is possible to assume that J « cV x B/4tt. 
However, because our numerical technique is relativist ic 
(see $3]), we cannot neglect the displacement current in 

Equation HI since it is used to evolve E. We thus choose 
to integrate Equation |4] neglecting the last two terms on 
the right hand side. Even though this approximation is 
not expected to affect the MHD-scale dynamics of the 
plasma, the neglected terms can still formally be larger 
than the displacement current dE/dt, which is not 
accurately evolved regardless of the vo <^ c condition. 
The effect of this can be seen in Equation [TJ where 
the last two terms of the right hand side represent the 
appearance of extra electric charges in the rotating 
box. In particular, the term proportional to vo can be 

much larger than V • E if \E\/\B\ <C \v \/c. This is 
expected in the case of the small-box approximation, 
where the typical magnitude of the turbulence velocity 
(~ c\E\/\B\) is much smaller than \vq\. As a first 
approach to this problem, we will neglect the existence 
of theses extra charges, assuming that they do not 
affect the plasma microphysics. Thus, our simulations 
will solve the standard Maxwell's equation, with the 
additional forces due to gravity acting on each particle 
individually. 

In the rotating frame, the particles will experience 
Lorentz forces, plus the Coriolis and tidal forces; in the 
case of a Keplerian disk, these are given by the well 
known expressions: 



2. BASIC EQUATIONS 

We carry out our study in the local, small-box approx- 
imation, where the size of the simulation box is much 
smaller than its distance to the center of the disk, r*o. 
The box rotates with the disk at the local (Keplerian) 
orbital frequency ujo = ^(^o)? so the reference frame is 
non-inertial. In the rotating frame, Maxwell's equations 



— - = 6muor ) xx 
at 



2ujq x p- 



■q{E+-xB), 
c 



(5) 



where p and u are the particle's momentum and velocity, 
m and q are its mass and charge, and x corresponds to 
the radial coordinate. 

The expressions for the Coriolis and tidal forces pre- 
sented in Equation [5] are valid in the "cold" limit (\u\ <C 
\vo |). Even though in our simulations the particles will 
reach relativistic velocities, the validity of the cold limit 
will still hold, but in a fluid sense. This means that, as 
long as the fluid velocity for each species satisfies \u\ <C 
\vq\, the fluid dynamics will be well described by the cold 
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limit expression (Equation [5]). Also, since in the MRI 
turbulence \u\ can be similar to va (= B / ^/Aitpi, where 
the subscript i stands for ions), our non-relativistic, cold 
limit will be strictly valid when va <C c. 

3. SIMULATION SETUP 

Our simulations are performed usi ng; the e l ectro- 
magnetic P IC code TRISTAN-MP (jBunemanl fl99l 
Spitkovsky 2005) in one and two dimensions. In 2D, the 
simulation box consists of a rectangle in the x — z plane, 
where x corresponds to the radial coordinate and z rep- 
resents the vertical direction of the disk. The azimuthal 
direction (into the simulation plane) is given by the y 
axis. The shearing velocity is v = —xsy, where s is 
the shearing parameter (= 3cjo/2, in the Keplerian case). 

In standard MHD simulations, shearing periodic bound- 
ary conditions are used along the radial (x) direction (see, 
e.g., Hawley et al. 1995). In that case, Galilean trans- 
formations of the MHD quantities at the boundaries are 
made to compensate for their initial velocity difference 
\Av\ = sL x . These shearing periodic boundary condi- 
tions can not be self-consistently implemented in PIC 
simulations. This is because, under a relativistic change 
of reference frame, the transformation of J cannot be 
obtained by only transforming the velocity of particles. 
This can introduce charge conservation problems at the 
box boundaries, even if \Av\ = sL x <C c. We avoid this 
difficulty by implementing shearing coordinates, in which 
the grid moves with the shearing velocity v = —xsy. In 
this new frame, the net plasma velocity at the bound- 
aries cancels out, allowing the use of periodic boundary 
conditions. In the shearing coordinate system, the 2D 
version of Maxwell's equations get modified by the pres- 
ence of extra terms in Faraday and Ampere's equations. 
The new equations read 



dB(r,t) 
dt 

dE(r,t) 
dt 



-V x E(f,t) - sB x (r,t)y and (6) 
V x B(r,t) - 4ttJ- sE x (r,t)y. (7) 



(see Appendix lAl for the derivation of equations [6] and [7]). 

Apart from the modifications to Maxwell's equations, 
forces on the particles will also transform in the shearing 
coordinate system: 



dp 
~dt 



2uj p y x - -ujqPxV + q(E 1-xB), (8) 



where p is the particle momentum. We can see that the 
combination of tidal and coriolis forces are substantially 
modified, with no dependence on the x coordinate in 
the shearing frame Q (see Appendix [A] for the derivation 
of Equation [8]) . Thus, our initial set up will consist of a 
periodic box where, apart from Lorentz forces, particles 

1 Notice that the combined expression for the coriolis and tidal 
forces in Equation [8] and the modified expression for the induc- 
tion law given by Equation [6] are equiv alent to the 2D versions 
of equations 14 and 15 of Johnson et al. ( 2008). These equations 
correspond to the MHD momentum evolution and induction equa- 
tions, expressed in terms of A v = v - v or ^, with v being the total 
fluid velocity and v or ^ = —3xouo/2y. 



are pushed by the forces corresponding to the first two 
terms in the right hand side of Equation [8j and where 
the fields are evolved according to Equations [6] and I7FI 

Our simulations are defined by a series of parameters, 
which set both the physical conditions and numerical 
resolution of the runs. The physical parameters are the 
ion to electron mass ratio m^/m e , the initial magnetic 
field direction and strength, the orbital frequency cjo, 
the initial ion and electron pressures pj, and the x and 
z sizes of the box (L x and L z ). The initial magnetic 
field along z (B Zy o) is quantified using the corresponding 

Alfven velocity of the plasma, v z A = 
The box size is normalized by Ao = 2irv A /ojo (roughly 
the wavelength of the fastest growing MRI mode in 
the MHD limit), and the orbital frequency is expressed 
in terms of the initial ion cyclotron frequency uj z ci 
(= \e\B Z fi/rriic), so our free parameter is the plasma 
magnetization uo^/ujq. Finally, the initial pressure of 
the particles is expressed in terms of their initial beta 
parameter, f3 z (= 8irpj/B^ Q ). 

In order to understand our choice of parameters, it is 
useful to know how they affect the total computing time, 
T comp , of the runs. The computational cost of a simu- 
lation is proportional to N ppc x N ts x N gp , where N ppc 
is the number of particles per cell, while N ts and N gp 
are the number of times steps and of grid points of the 
runs, respectively. Thus, it is possible to show that the 
computing time necessary to run a 2D simulation for a 
given number of orbital periods, Pq (= 2tt/ujo), scales as 



L comp ^ 



(miK) 3 / 2 ( c /^ i0 )Ky Wo ) 3 x 

\(c/oj p , e )/A x )\L/\ ) 2 N ppc (A x /(A t c)) , (9) 



where c/uo p ^ e is the electron inertial length, and and 
A t represent the grid spacing and simulation time step, 
respectively. Equation [9] shows that the computing time 
increases for large values of the mass ratio (rrii/m e ) and 
magnetization (uj z ^/cjo)? an d for small values of the ini- 
tial Alfven velocity (v z AQ /c). In addition, there is the 
increase in computing time due to spatial resolution 
(c/u) Pi e/A x ), box size (L/Ao), and particle resolution 
(Nppc)- Thus, in general, our simulations will use rather 
low values for the ion to electron mass ratio, m^/m e , and 
magnetization, cj^/gjo? and high values of v A0 /c. The 
low mass ratios and magnetizations used in our runs will 
be significantly far from realistic values. We will assess 
the dependence of our results on these parameters. 

4. ID SIMULATIONS: DISPERSION RELATION ANALYSIS 

In this Section we study the linear behavior of the col- 
lisionless MRI, and compare our results with previous 
analytical studies. We use ID simulations, where the x 
(radial) dimension is reduced to a few cells; this way only 

2 Since the modification to Faraday's (Ampere's) equation only 
affects the evolution of B y (E y ) with an extra term that depends 
on B x (E x ), we integrate B y (E y ) using simple time and space 
interpolations of B x (E x ). This way, after this modifications are 
implemented, the numerical algorithm used by TRISTAN-MP con- 
tinues to be second order accurate in time and space. 
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TABLE 1 

Parameters for the different sets of ID runs 



Runs ft? By i0 /B z , v a,o/ c ^c^/^o rni/m e 

01 0.05 1/20 11 10 

02 0.05 1/20 -11 10 

03 0.05 1/20 33 10 

04 0.05 1/20 -33 10 

05 0.05 1/20 110 10 

06 0.05 1/20 -110 10 

07 1 1/20 33 10 

08 1 1/20 -33 10 

09 10 1/20 33 10 

010 10 1/20 -33 10 

011 0.05 1/20 33 1 

012 0.05 1/20 33 5 

013 0.05 1/20 33 20 

014 0.05 1/60 33 10 

015 0.05 1 1/20 33 10 

016 1 1 1/20 33 10 

017 10 1 1/20 33 10 

018 1 1/20 220 10 

019 1 1/20 -220 10 

020 1 1/5 33 10 

021 1 1 33 10 

022 1 1/60 33 10 

Note. — We list the beta parameter of ions and electrons /3| 

(where j stands for ions or electrons), the ratio between the mean 
y (azimuthal) and z fields B y ,o/B z ,o, the initial Alfven velocity 
V A o/ c ' ^ ne plasma magnetization denned as the ratio of the initial 
ion cyclotron frequency and the disk rotation frequency ou^/ouq, 
and the ion to electron mass ratio m^/m e (/3|, v z A , and uj z i are 

calculated using only B Zj o). The space and particle resolutions 
in all of our ID simulations are given by c/w P) e/Aj; = 10 and 
Nppc — 15. 

wave vectors, fc, pointing along the z— axis are resolved. 
The box length along z, L z , is varied in such a way that 
only one single mode (with |fc| = 2tt/L z ) can grow. The 
mode is seeded by means of an initial plasma velocity 
^seed = (v^ /20) sm(27Tz/ L z )x^ which, by itself, would 
induce an Alfven wave of linear amplitude in the plasma 
(\SB\/B Zy o ~ 1/20). By measuring the growth rates 
in each case, we calculate the MRI dispersion relation, 
which then we compare with previous analytical results. 
The simulation parameters for each dispersion relation 
studied are specified in Table [TJ We ex plore bo th the 
weak field regime (jKrolik fc Zweibelll2QQ€t lFerrardl2QQ^. 
and also the high magnetization limit (jQuataert et al.l 
2002). Although our main focus is the highly magnetized 
case (lj^Jljo ^> 1), the study of the low magnetization 
limit will help us understand the parameter regime that 
minimizes the effect of weak fields, while optimizing the 
use of computer time. 

4.1. Low Magnetization Regime 

The low magnetization regime has been explored ana - 
lytically both in the cold limit (IKrolik fc Zweib el 2006), 
and in the finite temperature case (|Ferraroll2007l ). We 
will investigate the effect of uo^/uoq in the cold case first. 
We measure the MRI dispersion relation for /3| = 0.05 
using the magnetizations uj^/uoq = 11, -11, 33, -33, 110, 
and -110 (corresponding to simulations 01-06 in Table 
[T|) , which are presented in Figure [TJ The red lines cor- 
respond to Iuj^JujoI = 11 with the solid (dashed) line 
depicting the uj^Jujq > (< 0) case. Although nei- 
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Fig. 1. — Dispersion relations for ID simulations 01-06 from Ta- 
ble [TJ which use the same parameters ((3 Z = 0.05, v z A Q /c = 1/20, 

rrii/m e = 10, and B Vi o = 0), except for the initial magnetization 
of the plasma, which is denned by the ratio of the initial ion cy- 
clotron frequency to the disk orbital frequency uo z i /ujQ. The red, 
green, and black lines are for \uo z i /uJo \ = 11, 33, and 110; solid and 
dashed lines show lj^/ljq > and < 0, respectively. The results 
converge at \uj* i /uJo\ = 33 to a dispersion relation reasonably in 
agree ment with the analytical MHD prediction (Balbus & Hawlev 
1991), shown with the blue line. The growth rate uo is normal- 
ized in terms of the orbital frequency cuo, and the wavenumber k 
is normalized in terms of ko (= ujq/v z a q ). These numerical results 
are co nsistent wit h the a nalytical dispersion relation calculation of 
IKrolik & ZweSEel ([2006). 



ther the wave number of the fastest growing mode nor 
the corresponding growth rate varies between these two 
cases, the range of unstable wave numbers extends to 
larger values when lj^/ujq = —11. The green and black 
lines show results for |cj^/o;o| = 33 and 110, respec- 
tively. There is practically no difference between these 
two magnetizations, and the sign of uo^/uoq no longer 
plays a role. This shows that, when \uj^ i /ojQ\ = 33, our 
simulations have already converged to a dispersion rela- 
tion reas onably in agreement with the analytical MHD 
result (Balbus' & Hawlev 1991). The way the dispersion 
relation depends on the sign and magnitude of uj*Jujq 
shows that at low magnetization the coupling between 
the particles' gyromotion and their epicycle motion can 
modify the MRI dynamics significantly. These results 
are consistent with the MRI dispersion relations at low 
magnetizatio n (zero tempe r ature ) presented in Figures 1 
and 2 of IKrolik fc Zweibel (|2006h . 

We have also studied the effect of finite particle temper- 
atures (and, thus, of finite Larmor radius) by re-running 
simulations 03 and 04 using /3J =1 (runs 07 and 08) 
and 10 (runs 09 and O10). The results are presented 
in Figure El The black lines show the 'cold' (/?J = 0.05) 
case of Figure [TJ and the green and red lines show the 
/3f e = 1 and 10 results. The solid and dashed lines show 
oj*Ju)q > and < 0, respectively. Figure [2] shows that 
for larger initial plasma pressure, the range of unstable 
MRI modes shifts to larger wavelengths. The maximum 
value of the growth rate also increases for larger initial 
pressure. In addition, no substantial difference is ob- 
served between the uj z ci juj§ > and < cases. This 
result can be c ompared with the analytical treatment of 
iFerrarol (|2007| ). where the effect of finite Larmor radii 
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Fig. 2. — The black, green, and red lines show dispersion rela- 
tions for ID simulations with the same initial plasma magnetization 
\uj* i /(jJo \ = 33, but with different plasma betas: /3| = 0.05 (03 and 
04), 1 (07 and 08), and 10 (09 and O10); solid and dashed lines 
show uj z c i /ouo > and < 0, respectively. As /3| increases, the un- 
stable MRI modes shift to larger wavelengths, with the maximum 
growth rate increasing. This is a consequence of the ion Larmor 
radius increasing with increasing /3| . There is no substantial differ- 
ence between uj z c i /ouo > and < 0. The blue lines show runs like 
07 and 08, but with 6 times larger magnetization (\uj* i /a;o| = 220; 
runs 018 and 019). The migration to longer wavelengths seen in 
runs 07 and 08 is reduced, and the dispersion relations approach 
the (3j = 0.05 cases. This is consistent with the fact that relatively 
large values of f3j do not produce significant finite Larmor radius 
(FLR) effects if the magnetization is large enough. 

(FLR) of the ions was included. Ferraro's analytical pre- 
diction is consistent with our numerical calculation only 
for ojci/ojQ > (see his Figure 1). For uj^/ojq < 0, his 
prediction is that larger temperatures increase the upper 
limit of the unstable wave numbers and reduce the maxi- 
mum growth rate. We do not find this dependence on the 
sign of B z ^. This may be because our simulations use 
electrons in energy equipartition with the ions, in con- 
trast with the perfectly cold electrons used by iFerrarol 
(2007). In any case, if the dependence on /3| in Fig- 
ure [2] is caused by FLR effects of the kind predicted by 
IFerrarol (j2007l ). then increasing the plasma magnetiza- 
tion for fixed /3J should cause our results to approach the 
limit in which FLR effects are negligible. We tested this 
by taking runs 07 and 08 (with /3? e = 1) and increasing 
their magnetization by a factor of 6 to Icj^/^oI = 220 
(while keeping the same /3J). The corresponding disper- 
sion relations are shown by the blue lines in Figure O 
Increasing the magnetization indeed increases the range 
of unstable wave numbers, with the results approach- 
ing the cold plasma (zero FLR effects) results. Thus, for 
the magnetizations utilized in this paper, we expect FLR 
effects to be present for finite values for f3j. In the as- 
trophysical regimes of interest, however, FLR effects are 
expected to be negligible. 

4.2. High Magnetization Regime 

The linear behavior of the MRI has been studied analyt- 
ically in the high magnetization regime by ? using the 
kinetic MHD approach. In this approach, it is assumed 
that oJci/ojQ — ^ oo and also that FLR effects are unim- 
portant, i.e. the ion gyroradius is much smaller than 
Aq. One of the main differences with respect to stan- 
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Fig. 3. — Dispersion relation calculations for runs with B y ^ = 
B z ^o but different values of /3| = 0.05, 1, and 10 (black, green, 
and red lines corresponding to runs 015, 016, and 017, respec- 
tively). The growth of longer wavelength modes is favored at higher 
pressures, with the maximum growth rate for the /3| = 10 case be- 
ing significantly larger than for B yj o = (shown by the red lines 
in Figure |2 [). This is consistent with linear analytic calculations 
( Qua taert etalll^OOl ) . 

dard MHD is the increase in both the growth rate and 
the wavelength of the fastest growing mode, which hap- 
pens for large /3J in the presence of a significant toroidal 
magnetic field. We tested this result for simulations with 
/3J = 0.05, 1, and 10, with B yi0 /B z , =1 (runs 015, 016, 
and 017, respectively). The corresponding dispersion 
relations are shown in Figure [3j where the black, green, 
and red lines depict the cases with /?J = 0.05,1, and 
10. The tendency to favor the growth of longer wave- 
lengths and for larger maximum growth rates for larger 
/3j is clearly seen. Notice that the maximum growth rate 
for the /3j = 10, B y $ = B z $ case is significantly larger 
than for the analogous case with B y $ = (shown in 
Figure [2]). These dispersion relations are in reasonable 
agreement with their anal ytical counterpa r t show n in the 
right panel of Figure 4 of lOuataert et al.l (|2002[ ). 

4.3. Dependence on other Parameters 

Given the low mass ratios rrii/m e that we use, it is im- 
portant to check that our results are not affected by this 
parameter. Thus, we carried out simulations analogous 
to run 03 but using different values of rrii/m e . The dis- 
persion relations are shown in Figure [4] for rrii/m e = 1,5, 
10 and 20 (green, blue, black, and red lines, respectively). 
Figure [4] shows that rrii/m e does not play any significant 
role in the linear dispersion relation of the MRI, which 
does not change between rrii/m e = 1 and 20. Finally, in 
the same figure, a test of the effect of v z A /c is shown by 
the dashed-black line, which shows the dispersion rela- 
tion for v z A /c = 1/60 (run 014). There is no substantial 
difference relative to the v z A Q /c — 1/20 cases. 

4.4. MRI Saturation in ID 

Most of our non-linear MRI analysis will be done using 
2D simulations. In this section, however, we determine 
a saturation criterion for the magnetic amplification in 
ID. Although far from realistic, this will help us better 
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cases. 
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Fig. 5. — The Alfven velocities v A /c calculated with the k com- 
ponent of the magnetic field (k = x,y, and z) in three ID sim- 
ulations (022, O20, and 021). The runs have the same box-size 
L z = 1.25Ao (Ao = 2ttv a q/ljq) and the same parameters, except 
for their initial v z A /c. The black, green and red lines correspond 
to v AQ /c = 1/60, 1/5, and 1 (runs 07, O20, and 021 of Table [Q 
respectively). All three calculations have roughly the same satura- 
tion, with v x A rsj c and v v A ~ 10c. Similar results are obtained using 
a larger box (L z = 5Ao), and at different mi/m e and magnetiza- 
tions. 

understand and interpret the 2D results in ^5] Figure [5] 
shows the magnetic energy evolution for three ID sim- 
ulations with "box" sizes L z = 1.25Ao and with similar 
parameters except for their initial v A0 /c. The black, 
green and red lines correspond to v A0 /c = 1/60,1/5, 
and 1 (runs 022, O20, and 021 of Table [Q), and the 
solid and dashed lines show the evolution of their radial 
and azimuthal magnetic energies, B^/Bq and By/B^, 
respectively. The maximum amplification of B x in all 
three cases satisfies v A ~ c (where v A corresponds to the 
Alfven velocity calculated only with the x component of 



the magnetic field). After the saturation of the radial 
field, By continues to grow, but at a significantly lower 
rate. This result appears to be independent of the size 
of the box (it was also tested for L z = 5Ao), and other 
parameters like rrii/m ei u z ci ju§, and /3j. Note also that 
for v A /c = 1 the linear growth rate is reduced, leading 
to a suppressed exponential growth of B x . 
It is important to emphasize that, by assumption, our 
treatment of the MRI is valid only in the non-relativistic 
regime, i.e., when va <C c. This regime is the most 
interesting since va ~ c corresponds to a magnetic field 
energy close to the rest mass energy of the particles, 
which should be precluded by energy conservation, 
except very close to a black hole event horizon. Thus, 
this saturation criterion implies that, at least in ID, 
there is no mechanism stopping the growth of the field. 

In the next section we study how this ID evolution is 
modified by 2D effects. We also analyze the sources of 
angular momentum transport in detail, and study the in- 
terplay between the non-linear MRI turbulence and par- 
ticle heating. 

5. TWO-DIMENSIONAL SIMULATIONS 

Our 2D analysis is organized in four parts. First, ^5.11 
describes the overall non-linear behavior of the MRI tur- 
bulence, paying special attention to its saturation. £15.21 
briefly discusses how the non-linear evolution is modi- 
fied in the zero magnetic flux case. £15.31 analyses angular 
momentum transport, considering the contribution of an 
anisotropic pressure stress. Finally, in §5.41 we discuss 
particle heating and identify the different processes that 
contribute to it. 

Our analysis is based on a series of simulations listed in 
Table O The initial physical conditions of the runs are 
defined by: the beta of ions and electrons /3J, the mag- 
netic field along y, B y $, the field along £, B z $ (quantified 
via v A0 /c), the plasma magnetization lj^/ljo, and the 
ion to electron mass ratio, rrii/m e (v A ^ /c, and uj z c ^ 
are calculated only considering B z $). The remaining pa- 
rameters determine the numerical resolution of the runs. 
These are defined by: the box dimensions (L x x L z )/Xq 
(where, as before, Ao = 2ttv a /^o), and the space, time, 
and particle resolutions. The space and time resolutions 
are set by the number of grid points per electron skin 
depth, c/ujp }e /A x , while the particle resolution is defined 
by the number of particles per cell, N ppc . 

5.1. MRI Turbulence Evolution 
The non-linear MRI evolution is characterized by an ini- 
tial exponential growth of the field (until \B\/B Zy o ~ 5), 
followed by a significant decrease in the growth rate. 
This can be seen in Figure El which shows the evolution 
of the three magnetic energy components for simulations 
Tl (black) and T3 (red) of Table [2] quantified using their 
Alfven velocities (v A /c, v A /c, and v A /c, represented 
by solid, dashed, and dotted lines, respectively). In 
both cases, while v A /c continues to grow exponentially 
(although at a slower rate), the x and z components 
appear to enter a linear growth regime and saturate at 
amplitudes about one order of magnitude smaller than 
the azimuthal field. In both cases, the magnetic growth 
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TABLE 2 
Parameters of the 2D runs 



Runs 




By,o/ B z ,o 


V A,o/ C 




mi/me 




c/cu Pie /A x 


Nppc 


zero flux? 


rp-i 
1 1 


i 


U 


i /on 
1/zU 


1 1 


o 

z 


O X O 


1-7 

( 


Q 

o 


no 


rpr> 
1 Z 


i 


U 


i /on 
1/zU 


11 


r 




/I A 

4X4 


<-7 

( 


Q 

o 


no 


T3 


1 





1/60 


11 


2 


5x5 


7 


3 


no 


T4 


40 





1/20 


22 


2 


4x4 


7 


3 


no 


1 


1 


1 


1/20 


11 


2 


8x8 


7 


3 


no 


T6 


1 





1/20 


22 


2 


4x4 


7 


3 


no 


T7 


1 





1/120 


11 


2 


5x5 


7 


3 


no 


T8 


10 





1/20 


11 


2 


8x8 


7 


3 


no 


T9 


1 





1/20 


11 


2 


2x2 


10 


6 


no 


T10 


1 





1/20 


11 


2 


4x4 


7 


3 


no 


Til 


10 





1/20 


11 


2 


2x2 


14 


3 


no 


T12 


1 





1/60 


11 


2 


8x8 


7 


3 


yes 


T13 


1 





1/20 


11 


2 


8x8 


7 


3 


yes 


T14 


1 





1/20 


11 


2 


16 x 8 


7 


3 


yes 


Note. 


- A list 


of 2D simulations 


, denned by the initial beta parameter of ions and electrons /3J, the 


ratio between 


the mean y 



(azimuthal) and z (vertical) fields, B yj o/B Zj o, the initial Alfven velocity, v z A Q /c, the plasma magnetization, lj^/ljq (the ratio between 
the initial ion cyclotron frequency and the rotation frequency of the disk), and the ion to electron mass ratio, rrii/m e (the z superscript 
indicates that v A /c, and uo z i are defined by the z-component of B). The numerical resolution of the runs is defined by the box size 
L x /Xo and L z /Xo (where Ao = 2kv a /^o) 5 and the space and particle resolutions, which are determined by the number of grid points per 
electron skin depth c/uj p , e /A x and the number of particles per cell N ppc , respectively. 

plotted at t = 1.6, 3.6 and 14Pn- At t = I.6P0 the insta- 
bility is in the mildly non-linear regime (\B\/B Zj o ~ 1). 
At that point, magnetic amplification is dominated 
by the fastest growing MRI mode with a wavelength 
A ~ Ao- At t = 3.6Po, on the other hand, the field 
has been amplified to \B\/B Zy o ~ 60 and the dominant 
wavelength has grown to A ^ 4 An. This migration to 
longer wavelengths occurs along with the reconnection of 
magnetic field lines associated with the small wavelength 
modes as they become non-linear. Magnetic reconnec- 
tion takes place in thin current sheets, where B x and 
By switch sign. Along with the onset of reconnection, 
there is the appearance of loop- like structures (see, for 
instance, the one at (x,z) ~ (200, 200)c/oj Pie in Figures 
[TMZFO- The magnetic loops also correspond to regions 
of high plasma density, as can be seen from Figure [HI 
which shows the evolution of the plasma density, p/po, 
and B x /Bo for run Tl. The formation of loops is much 
clearer at t = 14Pq- At that stage, they appear as 
overdense regions (p/po ~ 10), in pressure equilibrium 
with the surrounding magnetic field. 

As By gradually decays, the magnetic energy corre- 
sponding to B x and B z stays rather constant. This 
energy is contained primarily in the magnetic loops, as 
can be seen from the field plots at t = 14Po in Figures 
[3-13 • At this point, the fields are in a quiescent state, 
with the loops experiencing almost no evolution (as 
can also be inferred from the smooth magnetic energy 
evolution after t > 13 in Figure [6]). This implies that, at 
late times, the growth of new linear MRI modes is dra- 
matically suppressed. This behavior can be explained by 
the migration of the growing modes to large wavelengths 
due to finite Larmor radius (FLR) effects at large /3J 
(see §4.ip . Indeed, for the relatively low magnetization 
LJci/ujQ = 11 of simulation Tl, the observed increase in 
temperature corresponds to /3J ~ 1000 (when only the 
initial P 2 ,o is considered), which makes the particle's 
gyroradii outside of the loops larger than Aq (after B y 
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Fig. 6. — Evolution of the three components of the magnetic 
field for 2D simulations Tl (v z AQ /c = 1/20; black lines) and T5 

(v A q/c = 1/60; red lines). Magnetic field values are expressed in 
terms of the Alfven velocities v A (= B^/^AiTp, with k = x,y, and 
z). The solid, dashed, and dotted lines represent the x, y, and z 
components, respectively. As in ID, the magnetic field saturates 
at v A z /c ~ 1 and v y A /c ~ 10. 



proceeds until v A « 2v A w c and v\ « 10c, which is 
qualitatively consistent with the ID saturation criterion 
found in §4.41 Indeed, there are only two important 
differences relative to the ID case. The first one is 
the growth of B Zl which can not occur in the ID case 
(V x E can not have a component along z, given that 
d x = d y = 0). The other difference is that, regardless 
of the initial v A0 , the exponential field growth stops 

at \B\/B Zj q ~ 5. This can be explained by the role of 
magnetic reconnection in dissipating the field energy in 
the non- linear regime of our 2D runs. Figure [3 shows the 
2D structure of the magnetic field for run Tl, with the 
three field components along with the magnetic energy 
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B x /B Zi0 B y /B Z|0 B Z /B Z|0 log(B 2 /B; i0 ) 




x/[cAV e ] X/[C/ Wp , e ] X/[C/ Wp ,.] X/[C/C p , e ] 

B X /B Z|0 B y /B Z|Q B Z /B Z|0 log(B 2 /B^ ) 



-60.22 60.22 -52.31 52.31 -64.16 64.16 -3.72 3.72 




200 400 600 200 400 600 200 400 600 200 400 600 

x/[c/c p , e ] x/[c/« p , e ] x/[c/« p , a ] x/[c/« p , e ] 

B x /B Zi0 B y /B Zi0 B Z /B Z|0 log(B 2 /B 2 , ) 



-1 10.19 110.19 -74.42 74.42 -106,45 106.45 -4.09 4.09 




x/[c/<y p , e ] x/[c/o; p , e ] x/[c/ Wp , e ] x/[c/ Wp ,J 



Fig. 7. — The three components of the magnetic field B^, £> y , and and its energy normalized in terms of the initial field B^o, for 
run Tl at three different times. The arrows in the \og(B 2 /B 2 ) plots show the projection of the magnetic field direction on the x — z 
plane. At t = I.6P0 (top row)the MRI is in the mildly non-linear regime, with the magnetic fluctuations dominated by the fastest growing 
linear MRI channel mode. At t = 3.6Po (middle row), the growth has migrated to longer wavelengths, with the dissipation of the short 
wavelength modes dominated by magnetic reconnection. Finally, at t = 14Po (bottom row), the turbulence is in a quiescent state, with no 
MRI modes and magnetic field concentrated in loops. At that point, growing MRI modes have wavelengths larger than the box size. 



has been significantly dissipated). Thus, FLR effects 
should significantly increase the wavelength of the 
unstable modes, presumably to length-scales larger 
than the box size. In a realistic astrophysical scenario, 
however, the MRI will not be suppressed by FLR effects. 
Indeed, uj^Jujo is typically many orders of magnitude 
larger than the value used in ou r simulations, which 
would make FLR effects negligible (Fer rarol[2QQ7| ). 

Also, in all of our simulations the quiescent state hap- 
pens after the MRI has reached the saturation condi- 
tion v x A « 2v z A « c and v\ « 10c. However, we expect 
the MRI to saturate before reaching this condition in a 
more realistic 3D problem. Indeed, we believe that the 
2D geometry of our simulations favors both the differ- 
ent evolution of B y (compared with B x and B z ) and the 
unrestricted field growth in the non-relativistic regime 
(va <^ c )- Loop formation makes reconnection a 2D 
phenomenon, which is favored if these structures are well 
resolved by the simulation. In our 2D runs, this is the 
case only for reconnection of field lines laying mainly on 
the x — z plane. The effect of the 2D geometry can be 



seen in Figure [8j which shows the magnetic energy evo- 
lution of 2D and 3D kinetic MHD versions of run Tl 
(us ing the same modif ied version of the ZEUS code as 
in iSharma et al.l 1200(f ). We see that the magnetic field 
energy evolves fundamentally differently in the 2D and 
3D runs. In the 2D case, the magnetic energy growth 
does not saturate, with the B y contribution dominating 
in the non-linear regime (like in our runs); while in the 

3D case the field saturates at \5B/B Z $\ ~ 10 with B y 
contributing nearly the same energy as the other two 
field components. 

5.2. The zero net flux case 

The MRI evolution presented above can be substantially 
modified if the net magnetic flux along z is zero. Figure 
[10] shows the density and B x of run T12, which is analo- 
gous to run T3 but with Bo = — s'm(x/L x )B Zy oz. In the 
zero net flux case the MRI initially grows faster in the low 
field regions. This is due to finite Larmor radius (FLR) 
effects, as explained in £ 14.11 When \Bq\ is small, the beta 
of the plasma is large, which increases the growth rate 
of the fastest growing mode (see Figure [2]). This effect 
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Fig. 8. — Comparison of magnetic energy evolutio n in 2D (black) 
and 3D (red) runs that use the kinetic MHD model of Sharma et al. 
(2006). These fluid simulations are analogous to simulation Tl, 
i.e., they use the same initial (3 and box size L XjZ /Xo (with the 
3D case having L y = L x , z )- The x (radial), y (azimuthal), and z 
(vertical) components of the magnetic energy are shown with solid, 
dashed, and dotted lines, respectively. In contrast to the 3D case, 
the 2D case does not saturate, qualitatively reproducing the lack 
of saturation in our 2D PIC simulations when va < c. Also, in the 
2D case the magnetic energy is dominated by the azimuthal field 
component (black-dashed line), which is also in agreement with our 
PIC results. This suggests that the 2D geometry of our runs plays 
a crucial role in precluding field saturation in the v A < c regime. 

also appears to be stronger when uj^Jujo < 0, which is 
consistent with the slightly different growth rates seen in 
Figure [2] for different signs of lj^/ljo. This FLR effect is 
not expected in realistic astrophysical settings. 

The MRI saturates at smaller amplitude with no net 
flux, compared with finite net flux. This can be seen in 
Figure [TTJ where the magnetic energy evolution of runs 
T12 (v z A /c = 1/60) and T13 [v z A ^/c = 1/20) are de- 
picted in black and red lines, respectively. With no net 
flux, the saturation is no longer characterized by a partic- 
ular value of va ~ c (as in the finite flux case; see Figure 
[6]). Instead, different values of v A /c saturate with sim- 
ilar amplification factors: B x /B z $ ~ 10, B y /B Zj o ~ 30, 
and B z /B z $ « 4. Also, the x-size of the box L x /Xo 
appears to play a role in the final saturation, as can be 
seen by comparing runs T13 and T14 (which are equal 
except for having L x /\q = 8 and 16, respectively). In- 
deed, T14 produces somewhat larger values of B\ than 
T13. This can be understood by noting that, for suffi- 
ciently large L x , regions of positive and negative initial 
B z will behave as spatially distinct regions, being able to 
reach saturation at amplitudes similar to the finite flux 
cases. Notice also that the earlier saturation suppresses 
the formation of strong channel flows and magnetic loops 
(see the density and B x configurations in the saturated 
state in Figures fTUl andHDjf). 

5.3. Pressure Anisotropies and Anisotropic Stresses 

The growth of a pressure anisotropy with respect to the 
local magnetic field can contribute to angular momen- 
tum transport via an anisotropic pressure stress, A xy j = 
—ApjB x B y /B 2 , where Apj = p±j— p\\j a nd j stands for 
ions and electrons (jQuataert et al.l 120021 : iSharma et al.l 
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Fig. 9. — The plasma density, p/po, and B x , for our fiducial 
2D simulation Tl at three different times. The arrows in the 
B x plot represent the magnetic field direction on the x — z plane. 
At t = I.6P0 (top row) the instability is in the mildly non-linear 
regime. At t = 3.6Po (middle row), the short wavelength modes 
have dissipated by magnetic reconnection and growth has migrated 
to longer wavelengths. Finally, at t = IAPq (bottom row), the tur- 
bulence has died away and growing MRI modes have wavelengths 
larger than the box size. In this quiescent state, both the magnetic 
field and plasma are concentrated in loops created by reconnection. 



2006). This pressure anisotropy is expected in regions 
where B is being amplified by the MRI, due to the adi- 
abatic invariance of the magnetic moment of the par- 
ticles, fjLj = p±j/pjB. The anisotropic stress may be 
comparable to the Maxwell stress, M xy = —B x B y j^rK in 
low-collisionality accretion disks, as was found by previ- 
ous non-linear studies of the MRI in collisionless plas- 
mas (Sh arma et al.l 120061 l2007f ). These studies used a 
fluid based approach, which did not evolve the pressure 
in an entirely self-consistent way. Instead, an approxi- 
mate model was used both to close the fluid equations, 
and to limit the growth of Apj . As mentioned in §T\ Apj 
is regulated by plasma microinst abilities (mirror, ion cy- 
clotron, electron whistler, etc) acting on scales compara- 
ble to the gyroradii of the different species. This provides 
a mechanism for pressure isotropization in the absence 
of Coulomb collisions. The effect of these instabilities on 
particles' velocities is kinetic in nature and can not be 
consistently captured by a fluid approach. Thus, to date, 
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Fig. 10. — The plasma density, p/po, and for 2D run T12 
at three different times, showing the overall MRI evolution in the 
case of zero net B z flux. The first, second, and third rows show the 
linear, non-linear, and post-saturation states at t = I.6P0, 2.4Po, 
and 5.6Po, respectively. With zero-net flux reconnection is much 
more vigor ous and leads to saturation of the MRI prior to va ~ c 
(see Figure [11} . Magnetic loops are also much less prominent. 
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Fig. 11. — The evolution of the three magnetic energy compo- 
nents for simulations T12, T13, and T14, with no vertical flux. The 
maximum magnetic field amplification (Bm a x/Bz,o) m runs T12 
(v AQ /c = 1/60; black line) and T13 (v AQ /c = 1/20; red line) ap- 
pear to be almost the same, despite their different initial Alfven ve- 
locity. This is fundamentally different from the finite B z flux cases, 
where the lower the initial v z A , the larger the amplification factor 

(so that at saturation v A ~ v z A ~ c, and v v A ~ 10c). The am- 
plification seems to increase in the case of run T14 (y z A Q /c = 1/20; 
green line), which has a L x /\q twice as large as the one in run T13. 
Thus, in zero net flux simulations, field amplification depends on 
box size (L X)2 /Ao), which is not the case in simulations with net 
vertical flux. 



Figure [12] shows the spatial distribution of the ion 
pressure anisotropy (Figure \l2h) and the correspond- 
ing anisotropic stress (Figure [T2b) at t — 1.6Pn,3.6Pn, 
and 14Po for run Tl. By comparing with Figure Hi, we 
see that at t = 1.6Pn the maximum anisotropy occurs 
in regions of large magnetic amplification, which conse- 
quently coincide with minima in /3\\^ (shown in Figure 
112b). The anisotropy Api/p\\^ is also well correlated 
with the Maxwell stress M xy (Figure fT2tZ) and, there- 
fore, with the anisotropic stress, A xy ^. At this early 
(linear) stage, the ion anisotropy satisfies Api/pu^ <C 1, 
thus plasma microinstabilities are not expected to pro- 
vide significant pressure isotropization (given that, as we 
will see below, microinstabilities isotropize the plasma 



their effect has been modeled by imposing a "hard wall" pressure efficiently when Api/p^ ~ l//3j* ., with q ~ 1, 
upper limit to \p±,j/p\\j — 1|, based on the assumption 
that Apj/p\\j will grow only until the relevant microin- 
stabilities reach their instability threshold. This crite- 
rion is motivated by solar wind observations ([Bale et al.l 
|2QQ9[ ). and theor etical and PIC s tudies of the relevant 
instabilities (e.g., Ga ry et aT] [l997). However, how these 
criteria apply given the simultaneous driving of the MHD 
turbulence by the MRI remains to be clarified. In this 
section, we describe the evolution of the anisotropy stress 
self-consistently using PIC simulations of the MRI, and 
quantify their contribution to transport in the disk. In 
£ 35.3.11 we provide a detailed 2D description of pressure 
anisotropics and their corresponding anisotropic stress 
A xy j. In £15.3.21 we analyze the dependence of Apj and 
A xy j on different simulation parameters using volume- 
averaged quantities. Finally, in ^5.3.31 we illustrate the 
growth of anisotropy-driven microinstabilities by identi- 
fying and analyzing the properties of the relevant small 
scale modes. 

5.3.1. Spatial distribution of the anisotropics 



and initially (3f = 1). The lack of ion isotropization can 
be seen in Figure [T5b . which shows that at early times 
the average magnetic moment of ions, remains very 
close to its initial value /i^o- 

At t = 3.6Po, on the other hand, the correlation 
between A xy ^ and M xy gets significantly suppressed. 
A XVi i is especially suppressed in regions of large M xy , 
which coincide with the regions of the lowest /3m }i . 
This suggest the presence of an efficient mechanism for 
pressure anisotropization at low f3\\j. We will quanti- 
tatively discuss (in §5.3.2p the most likely mechanism 
for ion pitch angle scattering in these low f3\\j regions. 
Magnetic reconnection is also expected to reduce Apj 
at this stage. Indeed, as we will see below, at this 
time reconnection contributes significantly to particle 
energization. Since this process is not expected to 
preserve /ij, pressure anisotropics must also be reduced 
due to reconnection. The non-conservation of fij can 
be seen explicitly in Figure [T3b . which shows that at 
late times /^//i^o ~ 60 on average (in accordance with 




Fig. 12. — 2D plots of the ion pressure anisotropy Api/p\\^ (Api = p±^ — P||,i) 5 the ion parallel beta the ion anisotropic stress 

A xy ,i/po, and the Maxwell stress M xy /po, at t = I.6P0, 3.6Po, and 14Po for run Tl (where po is the initial pressure in the plasma). Overall, 
the pressure anisotropy increases as the MRI grows. At late times, however, the pressure is roughly isotropic in the magnetic loops but 
anisotropic elsewhere. 
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Fig. 13. — 2D plots of the average magnetic moment of ions m 
(= p±,i/(piB)) at t = 1.6Po and 3.6Po for run Tl, normalized in 
terms of its initial value /i^o- At late times, the magnetic moment 
has increased significantly, consistent with the isotropization of the 
plasma pressure by small-scale kinetic instabilities (see Figures [15] 
andfl6l). 



B/B « 60 and f3 i: ± « 1). 

At t = 14Po there are no regions in the box where 
/3\\ t i <C 1, which suggests that the correlation between 
M xy and A xy ^ must be to some degree recovered. How- 
ever, this does not occur. A significant fraction of the 



magnetic field energy is contained in loops, but the pres- 
sure anisotropy within these loops appears to be almost 
zero. This is consistent with the fact that no signifi- 
cant magnetic amplification occurs in the loops, so they 
do not develop pressure anisotropics due to p conser- 
vation. Indeed, loops are a byproduct of the magnetic 
reconnection of the MRI- amplified field, so no significant 
growth of Api/p\\j is expected to occur in these regions. 
Thus, we see from Figures [T2fc and[T2t that, whereas the 
largest contribution to M xy comes from the inner part 
of the loops, the largest magnitudes of A xy ^ come from 
their outer parts. We also note that both M xy and A xy ^ 
can get large negative and positive values, implying that 
these quantities may produce a negative stress on aver- 
age. This is because, as can be seen from Figures [3 
and [7]/, the magnetic loops do not develop an anticorre- 
lation between B x and B y , as is the case with the MRI 
channel modes. As we will see below, this produces nega- 
tive volume- averaged M xy and A xy j in the late stages of 
some of our simulations. This is, we believe, an artifact 
of how our 2D simulations saturate. 



5.3.2. Volume- averaged anisotropies and stresses 
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Fig. 14. — The evolution of volume- averaged quantities for sim- 
ulations Tl and T4. The first row shows the different stresses 



M xy (red), R xy (R xy = R x 



R xy ,e\ green), and A xy (A xy 



A xy ^i + A xy ^ e ; black), normalized in terms of the perpendicular 
plasma pressure (p± = p± i +P_L, e )- The second row shows the 
parallel and perpendicular betas (in solid and dotted lines, respec- 
tively), for ions [black] and electrons [green]. These quantities are 
calculated dividing the volume averages of the particles' pressures 
and the magnetic field pressure. The third row shows the pressure 
anisotropy Api/pn^, of ions (black line; electrons follow a qualita- 
tively similar trend), along with the theoretically estimated thresh- 
olds for the ion-cyclotron (IC) instability (Ap^/cS green) and an 
empirical threshold obtained from solar wind ion anisotropy mea- 
surements (Api^sw), shown by the red line (Bale et al. 2009). Note 
that the plasma remains near the solar wind threshold (which is 
very similar to the threshold for excitation of mirror modes), par- 
ticularly in the higher beta simulation T4 (right column). 



In this section we quantitatively analyze the physics 
of pressure anisotropy evolution and the corresponding 
anisotropic stresses, using volume- averaged quantities. 
Figure [14] shows the time evolution of volume- averaged 
stresses, plasma betas, and ion pressure anisotropics for 
runs Tl and T4. The first column concentrates on run 
Tl, with panel [T4h depicting the Maxwell stress (M xy ] 
red), the Reynolds stress (R xy = R X yi + Rxy,e] green) 



and the anisotropic stress (A xy = A xy ,i + A xy ^\ black). 
These stresses are plotted as solid (dotted) lines when 
they are positive (negative), and are normalized in terms 
of the total perpendicular pressure p± = p± : i + p±, e - 
They are thus an estimat e of the contribution to th e 
effective a parameter of iShakura fc Sunvaevl (|l973). 
We see that in the exponential growth regime (from 



t = to t « 2Po) M xy and R xy almost coincide, in 
accordance with the expected linear MRI behavior. On 
the other hand, A xy appears to grow exponentially at a 
rate larger than that of M xy and R xy . This is because 
A xyj j/M xy ~ Apj/B 2 , so in the linear regime this ratio 
grows as Apj/B^ Q . Also, this ratio can be expressed as 
A xy j/M xy = Apj/3\\j/(2p\\j). Thus, since in the linear 
regime Apj/p\\j <C 1 and /3\\j = 1 (for run Tl), the 
ion and electron anisotropic stresses A xy j must satisfy 

A X y J <C M X y 

After t « 1.5Po, the growth of the anisotropic stress A xy 
becomes significantly slower than that of the Maxwell 
stress M xy . This can be understood in terms of the de- 
pendence of Apj/p\\j on /3\\j. Figure [14^ shows the time 
evolution of the ion pressure anisotropy Api/p^^ for run 
Tl (black line; the electrons follow qualitatively the same 
trend). When /3m^ < 0.3 (t < 3Po), ^Pi/p\\,i appears 
to be constrained by the condition for io n-cyclotron in- 
stabili ty growth (green line), as used by iSharma et al.l 
(12006). This threshold is given by 

^ < ^ (10) 



< 



which implies that A xy j/M xy = Apjf3\\j/(2p\\j) 

<C 1 if /3\\ t i <C 1. This explains the lack of 



0.18/?°;f 



correlation between A xy j and M xy in the low f3\\^ 
regions at t = 3.6Po, as seen in Figures [12^ and \V2h. 

At larger values of f3\\^ (t > 3Po), the pressure 
anisotropy appears to be limited by a condition less strin- 
gent than that of the ion-cyclotron instability. Remark- 
ably, the bound on the ion anisotropy in our simulations 
is consistent wi th the maximum anisotropy measured in 
the solar wind (|Bale et al.ll2009f ). shown by the red line. 
This limit is given by 



< 



0.77 



&pi 

P\\,i ~ 0%i + 0.016) 



0.76 * 



(ii) 



The bound in equation [TT] is very similar to the stability 
condition for the mirror instability, which is given by: 



Api < 1 



P\\,i 



(12) 



Thus both our simulations and the solar wind data sug- 
gest that the mirror instability plays the key role control- 
ling the pressure anisotropy when f3\\^ > 0.3. The signifi- 
cant role of the mirror instability will be confirmed below 
by directly identifying the presence of mirror modes em- 
bedded in the MRI turbulence. 

Figure [TH shows that the anisotropic stress begins 
to dominate for t > 6Po. In addition, at late times 
the pressure anisotropy Api/p\\ ti is significant larger 
than the limit provided by Equation [Til This can 
also bee seen in Figures [T2t and [T2]? '. which show the 
spatial distribution of Ap^/piu and /Su^ of run Tl at 
t = 14Po. Indeed, in most of the volume /3iu ~ 20, 
while Api/p\\i w 4 (see, for instance, the region centered 
at (x,z) = (250,500)c/co>p 5 e), which is significantly larger 
than what is expected from Equation [TU This lack of 
isotropization is probably due to the large value of the 
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ion Larmor radius Rl,% at the end of the simulations, 
which makes the typical wavelength of the mirror 
instability close to the MRI wavelength. Indeed, in 
regions of large anisotropy, Rl,i ~ 100c/ 'u p ^ while the 
wavelength of the dominant MRI mode is S00c/uj Pje (see 
Figure Wo)- Thus, given the similarity between mirror 
and MRI wavelengths, we do not expect the mirror 
modes to grow as effectively as they do in the regime 
where these scales are well separated^ 

The large beta behavior of the ion anisotropics can be 
seen in the second column of Figure [HJ which shows 
results for run T4, which initially has /3J =40. We see 
that initially A xy becomes comparable to M xy during 
the stage of exponential growth. Notice that A xy is 
similar to M xy only until the end of the exponential 
growth regime. After that, the significant decrease in 
makes M xy the dominant stress. We also see in 
Figure fT4b that M xy and A xy may acquire negative 
values. This indicates, as seen in §5.3. 1| that in the 
post-saturation state, stresses may be dominated by 
loop-like structures where B x and B y are not necessarily 
anti-correlated (as in the case of MRI modes). 
Figure [14] /' shows that the maximum of Api/p\\^ is 
determined by Equation [IT] almost all the time (except 
at the end of the run, as in run Tl). This indicates 
that using a large initial /3^j ensures that isotropization 

will be dominated by the mirror instability even at the 
end of the exponential growth, with the ion-cyclotron 
instability playing a less important role. This large 
beta behavior also happens for /3n = 10, as in runs 
T8 and 11. An interesting question is whether using 
/3J significantly larger than 40 would make A xy ^ the 
dominant viscous stress in the MRI-driven turbulence, 
because the precise value of A xy ^ depends on the exact 
dependence of Api/p^^ on /3||^. Testing this possibility 
requires using values of lj^Jljq significantly larger than 
the ones used here (in order to keep the ion Larmor 
radii much smaller than the dominant MRI wavelength, 
Ao). This possibility will be investigated in a future work. 

We also studied the dependence of Apj/p\\j and the 
plasma stresses on other simulation parameters. We 
tested the dependence on: mi/m e (with run T2; using 
Trii/m e = 5), different values of v z A /c (with run T3 and 
T7; using v\ /c = 1/60 and 1/120, respectively), initial 
finite azimuthal flux with B Vj q = B Zj o (run T5), and a 
larger magnetization uj z ci juj§ = 22 (run T6). All these 
simulations (which have /3m j = 1 as in run Tl) have 
the same behavior as run Tl, showing that the physics 
of pressure isotropization obtained in our simulations is 
reasonably well converged, with /3\\j being the only rel- 
evant parameter. 

Finally, we have also tested numerical convergence using 

3 Since the ratio of the MRI growth rate and the ion cyclotron 
frequency in our simulations is ujq/uj c ^ ~ 0.1 — 0.01, the mirror 
and ion-cyclotron time scales are much closer to the MRI time 
scale than in reality (where uoq/uj c ^ ~ 10 -7 ). However, our lin- 
ear calculations show that the pressure anisotropy thresholds for 
a much lar ger growth rate are only a factor of ~ 2 larger than 
in Equation [TT] and can only partially explain the large departure 
from the threshold seen in Figure [T4fe . 



run T9 and T10. Run T9 tested box size dependence by 
using L x = L z = 4Ao (half the values used in Tl), and 
run T10 tested space, time, and particle resolution by 
using c/u Pye /A x = 10 and N ppc = 6. No difference with 
respect to the results of run Tl were found. 

5.3.3. Mirror mode analysis 

As seen in §5.3.2| the physics of pressure isotropization 
at /3||j > 0.3 appears to be well described by the mirror 
instability, both for ions and (the large mass) electrons. 
In this section, we analyze the structure of the small scale 
modes in the plasma, and check whether they satisfy mir- 
ror mode properties. Figure [T5l shows the case of run T8, 
which is analogous to Tl but with /3| = 10 (the larger /3J 
makes the effects of the mirror mode more prominent, as 
can be seen in Figure Elf). Plot[T5b shows \og(B 2 / B z $) 
at t = 2Po. It is clear visually that, in regions of ampli- 
fied magnetic field, small scale fluctuations arise. The 
length scale of these modes is about 2Rl^ (where Rl,% 
is the average Larmor radius of the ions), as can be seen 
from plots [T5b-[T5b. These plots depict \og(B 2 / B Zi q) and 
SBj/ < B > in a zoomed region marked by the small 
rectangle centered at x/Rl^ = 100 and z/Rl^ = 99 in 
plot [T5h (where SBk = B^— < Bk > is a measure of the 
magnetic field fluctuations along the k— axis, and <> 
represents volume average within the zoomed-in region) . 
The typical wavelength of ~ 2Rl^ of the mirror modes 
is also observed in simulation T4, where the initial ion 
temperature is four times that of simulation T8. Also, we 
confirmed the numerical convergence of this scaling using 
run Til, which has a spatial resolution of c/cu Pje = 14 
(twice the one of run T8). 

One of the properties of the mirror instability is the 
anti-correlation between B 2 and plasma density p. This 
anti-correlation is present in Figures [T6h andfTBb. which 
depict (B 2 - < B 2 >)/ < B 2 > and (p- < p >)/ < p > 
in the same zoomed-in region shown in Figure [TKl Also, 
mirror modes satisfy SB JL k. This can be checked by 
computing the components of SB parallel and perpen- 
dicular to fc, shown in plots [T6b and [TBH respectively 
(where k and x are estimated to form an angle of 120°). 
The amplitude of SB± is significantly larger than that of 
SB\\, consistent with the mirror mode polarization. Fi- 
nally, we can compare the projection of SB on the x — z 
plane (shown by the SB X and SB Z plots in Figures [T5b 

and fT5b ) with the analogous projection of B (shown by 
black arrows in Figure [T5b ) . The magnetic fluctuations 
in the x — z plane roughly satisfy SB\\B, another polar- 
ization property of the mirror modes. 

5.4. The Energy Distribution of Particles 

In this section we explore the signatures of the differ- 
ent heating processes in the energy spectra of the par- 
ticles by analyzing the case of run Tl. Figures [3 and 
[9] show that at t = 3.6Po efficient migration into longer 
wavelength MRI modes is happening, with a correspond- 
ingly significant rate of magnetic field energy dissipation 
through reconnection. Figure [TTb shows that the average 
ion and electron spectra at that moment are composed of 
a thermal distribution, plus a power-law tail with spec- 
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Fig. 15. — Magnetic energy and magnetic field fluctuations 
for run T8 highlighting the appearance of mirror modes at t = 
2Po, driven by plasma pressure anisotropics. Panel a shows 
log(£> 2 IB 2 Z ) and panel b shows the same quantity but in a zoomed 
region centered at x/Rlj = 100 and z/Rlj = 99 (marked by 
the yellow-dotted rectangle in panel a). Panels c, d, and e depict 
SBk/ < B >, with k = x,y, and z. Rl,i is the average ion Larmor 
radius and <> represents volume average within the zoomed-in 
region. The x and z axes are normalized by Rl,i because this is 
the characteristic wavelength of the fastest growing mirror modes. 

tral index of ~ 1.50 (black and green correspond to ions 
and electrons, respectively). On the other hand, Fig- 
ure [TTb shows the energy spectra at t — 14Pq? which 

4 Notic e that this is the sam e recon nection-driven spectral index 
found by Sironi & Spitkovskv (2011) in their PIC simulations of 
relativistic striped wind shocks. This suggests a possible connec- 
tion between relativistic and non-relativistic reconnection-driven 
non-thermal particle spectra. 




Fig. 16. — The fluctuations in magnetic ene rgy and plasma den- 
sity for the same region depicted in Figure ^] are shown in panels a 
and 6, respectively. The anti-correlation between these two quan- 
tities is a signature of the mirror instability. Panels c and d show 
magnetic field fluctuations that are perpendicular and parallel to 
the dominant mirror wavevector k, respectively. The components 
parallel to k have an amplitude smaller than the perpendicular 
ones, showing that these modes roughly satisfy the mirror polar- 
ization SB _L k. 

corresponds to the "quiescent" state, where no reconnec- 
tion happens whatsoever. We see that, at this stage, 
there is not a power-law tail and, instead, a high-energy 
bump appears. The high-enegy particles concentrate in 
the regions outside of the magnetic loops, where A xy j 
is relatively large, as can be seen from Figure [T2fc. This 
indicates that at later times particle energization occurs 
mainly via a viscous heating propo rtional to A X y^ (lik e 
the one suggested in equation 6 of iSharma eF al. 2007). 
The energy spectra between t = 3.6Po and t = 14Po 
are a combination of these two distributions. In a fu- 
ture study, we will explore significantly larger values of 
rrii I m e to shed light on the different heating mechanisms 
for ions and electrons in collisionless disks. 

6. SUMMARY AND DISCUSSION 

In this paper we have studied the magnetorotational 
instability (MRI) in a collisionless plasma using first- 
principles particle-in-cell (PIC) simulations. Our 
motivation is the application to low acc retion rate, ra- 
diatively inefficient accretion flows (e.g., iNaravan et al.l 
1998). These flows are expected to be present in systems 
accreting at less than a few percent of the Eddington 
rate. This includes the central black hole of our Galaxy 
(Sgr A*) and most nearby galaxies, as well as the 
low-hard state of X-ray binaries. 
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Fig. 17. — The particle spectra for run Tl at times t = 3.6Po 
(panel a) and t = IAPq (panel b). Black and green lines correspond 
to ions and electrons, respectively. While the early-time spectra are 
composed of a thermal (shown for reference in red-dotted line) and 
power-law distribution, the late-time spectra correspond to a two- 
temperature distribution for each species. The formation of the 
power-law distribution is clearly correlated with magnetic recon- 
nection, while the heating at later times appears to be dominated 
by viscous heating due to the anisotropic stress. 



In the first part of the paper (SI), we studied the linear 
dynamics of the MRI using ID simulations in which only 
wavevectors along the z-direction - the rotation axis - 
are resolved. We focused in particular on understanding 
the effects of the plasma magnetization uj*Jujq (the 
ratio of the initial ion cyclotron frequency to the disk 
rotation frequency) and finite particle temperature 
on the dispersion relation of the MRI. Our results 
were in reasonable agreement wit h previous ana l ytical 
studies of the collisionless MRI (iQuataert et al.l 120021 : 
iKrolik fc ZweiT^IM^fFerrard 

In the second part of the paper, we studied the 
multidimensional aspects of the MRI using local 2D 
(axisymmetric) simulations. Given the computational 
effort involved in PIC simulations, we defer fully 3D 
simulations to future work. In our 2D calculations, 
we focused on the evolution of the plasma pressure, 
magnetic energy, and plasma stresses, and assessed their 



dependence on the different physical and simulation 
parameters. We found that the overall evolution of 
the MRI in axisymmetric simulations is qualitatively 
similar to that found in previous MHD simulations and 
simulations that inc luded fluid models of kinetic effects 
(Sh arma et aLl 12006) . In particular, in the PIC simu- 
lations with net magnetic flux, the MRI only saturates 
when the Alfven speed becomes comparable to the speed 
of light (irrespective of the initial value of the Alfven 
speed). This is consistent with the absence of saturation 
of MRI channel models in analogous axisymmetric fluid 
simulations (see Figure [8]). By contrast, for simulations 
with no net magnetic flux, the MRI saturates at a lower 
amplitude before va ~ c (see Figure [TT]). 

By adiabatic invariance of the magnetic moment of 
particles, the amplification of the magnetic field by 
the MRI in a collisionless plasma produces pressure 
anisotropics, Apj = p±j — p\\j, for each species j 
(where the directions are measured relative to the local 
magnetic field). These anisotropics can be important 
for angular momentum transport and particle heating, 
since they give rise to an anisotropic pressure stress, 
A xy j = —B x B v Ap j/B 2 , that may be comparable to 
the m agnetic stress ( Quataert et al.l [200 2t iSharma et al.l 
l2006f ). This is effectively a macroscopic viscosity in 
a collisionless plasma. We investigated this physics 
in a self-consistent way by studying the evolution of 
the plasma pressure during the linear and non-linear 
phases of the MRI. We found that the importance of 
the anisotropic stress depends on the instantaneous 
plasma /?. For small ft < 10, the ion anisotropic 
stress, A xy ^, is smaller than the Maxwell stress M xy . 
However, A xy ^ may be comparable or even surpass M xy 
when ft > 10. This regime is difficult to achieve in 

our simulations due to the continuous build-up of 5, 
which grows until v x A and v z A ~ c, and v y A ~ 10c. This 
difficulty could in principle be overcome by using large 
initial $j. However, this would require using values of 
the magnetization parameter uj z ci juj§ significantly larger 
than those used here (so that the MRI and microinsta- 
bilities length scales are properly separated), which is 
computationally difficult (see Equation [9]) . We expect 
that the saturation of both the magnetic field and the 
pressure anisotropy will be fundamentally different in 3D 
simulations. In that case, magnetic field amplification 
is expected to stop when va <C c, which would allow 
the existence of a turbulent state in which ftu > 10 
(see, e.g., the 2D vs. 3D fluid simulations in Figure 
[Sj). Our results suggest that transport and heating due 
to the anisotropic stress will be important in this regime. 

Absent any mechanism of temperature isotropization- 
which is provided by collisions in the fluid limit - the 
pressure anisotropy would continue to grow throughout 
the linear phase of the MRI in a low-collisionality 
plasma. In our calculations, we find direct evidence 
for collisionless processes that provide temperature 
isotropization in the absence of Coulomb collisions. In 
particular, the p± > p\\ state that is generically created 
by the MRI is unstable to the excitation of mirror 
modes (see Figures fT5l and [T6|). Moreover, once mirror 
excitation commences the volume- averaged pressure 
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anisotropy remains near the threshold for the onset 
of the mirror instability, particularly in our higher 
/3j simulations (as shown in Figures [Ml and Figures 
rB|f). This is consistent with the evolution of the 
pressu re anisotropy ob served in the near-Earth solar 
wind (jBale et al.l l2009f ). Our results, analogous theo- 
retica l work in the solar wind context (jHellinger et al.l 
2006), and the solar wind measurements all strongly 
suggest that temperature isotropization in moderate j3 
low-collisionality plasmas is dominated by (reasonably) 
well understood velocity space instabilities (in partic- 
ular, the mirror, firehose, and ion cyclotron instabilities). 

In future work, we will study the heating of particles 
in MRI turbulence in detail, paying particular attention 
to how the assumed electron to ion mass ratio m e /rrii 
gives rise to different electron and ion heating physics. 
In our preliminary analysis in the present paper, we 
found that particle heating occurs via two dominant 
mechanisms: one is reconnect ion and the other is the 
viscous heating prod uced by the anisotropic stress A xy j 
(jSharma et al.l 120071 ). These different heating processes 
are expected to imprint different signatures in the energy 
spectra of particles. Perhaps most interestingly, we 
find that magnetic reconnection produces a distinctive 
power-law distribution function (with spectral index 
~ 1.5) as the MRI becomes nonlinear, but before the 
turbulence has saturated and died away (see Figure 
[T7b). This spectral index is similar to the one found 
by iSironi fc Spitkovskvl ([201 lh in their PIC study of 
reconnection forced by relativistic shocks in the striped 
wind of pulsars. Whether this non-thermal energization 
is due to the reconnection electric field parallel to the 
curre nt sheet (as in the case of ISironi fc Spitkovskvl 
l2011f ) , or by a Fe rmi- like process du e to the formation of 
magnetic loops ([Drake et al.l 1200(f ). will be clarified in 
detail in a future work. The importance of reconnection 
diminishes once the MRI saturates. At that point, the 
particle energization appears to be dominated by the 
anisotropic stress A xy j, which is largest outside of the 
magnetic loops. This process produces a high-energy 
bump in the energy spectra of particles (see Figure fTTb). 
with the most energetic particles concentrating outside 
of the magnetic loops. In our calculations we do not 
find unambiguous evidence for heating via the turbulent 
cascade of Alfven and slow waves that is expected 
to be set up via the nonlin ear saturation of the MRI 
([Quataert fc Gruzinovlll999[ ). This could be a limitation 
of our 2D simulations which do not develop sustained 
turbulence; in addition, such a cascade may be diffi- 
cult to resolve given the modest dynamic range in our 



simulations between the inner and outer turbulent scales. 

The PIC simulations presented in this paper have 
numerous advantages relative to fluid calculations for 
studying the physics of the MRI in low-collisionality plas- 
mas. In particular, our PIC simulations represent the 
first self-consistent study of magnetic field amplification 
and saturation, particle heating, and pressure anisotropy 
evolution in MRI-driven turbulence. There are, however, 
also drawbacks associated with PIC simulations. In 
addition to being very computationally demanding, the 
need to resolve the electron skin depth implies that 
there is always an unphysically small dynamic range 
between the initial ion cyclotron frequency uo c ^ (where 
i stands for 'ions') and the disk rotation frequency 
uoq\ our simulations initially have uj c ^/ojq ~ 10 — 100 
instead of u Cy i/uJo ~ 10 7 in real systems. This ratio, 
however, increases to u Cy i/uJo ~ 100 — 1000 in the 
nonlinear regime, considering both the magnetic field 
amplification and the growth of the mean Lorentz factor 
of the particles. We find no strong dependence of our 
results on the initial value of uj Cy i/uo, suggesting that 
we are adequately in the "MHD" regime, but this will 
need to be carefully studied in 3D simulations as well, 
where the computational requirements are even more 
demanding. 

Another limitation of the current calculations is that our 
results are strictly valid only in the limit <C c (where 
is the bulk rotation velocity of the plasma). In partic- 
ular, the PIC equations evolved here neglect additional 
terms in Maxwell's equations that arise from being in a 
rotating reference frame (see Appendix [A]) ; these are self- 
consistently small in the limit of small rotation velocities. 
Our analysis does consistently allow for relativistic tem- 
peratures and thus it is possible to study quantitatively 
the limit of relativistic electrons, but non-relativistic ions 
that is of particular interest for radiatively inefficient ac- 
cretion flow models. 
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APPENDIX 
SHEARING COORDINATES 

We describe the plasma in the shearing frame, S' ', which is related to the usual cartesian frame S by the following 
coordinate transformation: 

x' — x, v' — r(y — vt), 

z' = z, and t' = T(t-vy/c 2 ), 

where the primed coordinates correspond to the frame S' , T = (1 — v 2 /c 2 ) -1 / 2 is the Lorentz factor due to shear, 
v = — sx is the ^/-component of the shear velocity, and s is the shear parameter (s = 3cjo/2 in the case of a Keplerian 
disk). In S' the plasma is initially at rest (with no shearing velocity), which allows us to replace shearing box periodic 
boundary condition by standard periodic boundary conditions. 

The transformation of the electric and magnetic fields from S to S' can be expressed as the standard Lorentz invariant 
transformation 

E' v = E y , E' ± = T(E ± + f x B ± ), 

B' y = B V1 and B' ± = T(B ± -fx E±). V ; 

As we will see below, the evolution of the fields in S' is given by modified versions of Maxwell's equations. The 
modifications arise from the dependence of v and r on x. Although we are interested in the small box limit, we derive 
the evolution of E' and B' assuming an arbitrarily large box, with the small box limit taken only at the end of the 
calculations. Since the module of v(x) can not exceed c, we will use a shear profile where this is satisfied regardless of 
the value of x. To do that, we will impose that the "local" shear profile seen by an observer moving with the flow is 
viocai = —sxiocal, where xi oca i = at the observer's position. One can show that this condition implies a global shear 
in the box given that v/c = — arctanh(sx/c). With this space dependence of v(x), the x— derivatives of v and T will 
be given by 

dv/dx = -s/T 2 and dT/dx = -sTv/c 2 , (A3) 
which we will use in the derivation of the field dynamics described below. 

Modifications to Maxwell's equations 

We determine the changes to each component of Maxwell's equations one by one. Let us start with Faraday's equation. 

Faraday's equation: x component 

Given the transformations defined in Equations IA1I and IA21 we want to know how the equation 

^ = -(cVx£(f,^ (A4) 

is modified in the coordinate system S' . From Equations IA2I we get 

B' x {f,lf) = T(B x (f,t) - *E z {?,t)), E' z {fJ) = T(E z (f,t) - *B x {?,t)), and E' y ^,t') = E v (r,t) (A5) 
Then, from Equation IA1I it is possible to show that 

8B' x (f,e) _ r2 / dB x (r,t) _ v8EAr 1 t)_ +v dB x (r,t) v 2 0E z (r,t) 



and 



dt' V dt c dt dy c dy 

8E f y (f,t f ) _dE y (r,t) 



dz r dz 

8E' z (?,t f ) __ 2( v dE z (f,t) v 2 dB x (r,t) , dE z (f,t) v dB x (r,t) 



\c 2 dt c 3 dt dy c dy J 

(cV xE'(7,t')) x . (A7) 



dy' \c 2 dt c 3 dt dy c dy 

Combining these three expressions, one can show that 

dB' x {fJ) _ 
dt' 

So the x— component of Faraday's equation does not transform from S to S' . 
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Faraday's equation: y component 



Here the procedure is analogous to the case described above. However, since the y component incl udes partial 
derivatives with respect to x' ', the x' dependence of v and T will give rise to new terms. From Equations (jA2|) we get: 

B;(r<,0 = B y (f,t), E' x (?,t) = T{E z {r,t)-lB x {r,t)), and E' x (? J) = T(E x (f,t) + *B x {f,t)). 

(A8) 



Then 

i 



dB' y (r>,t>) ( dB y (r,t) dB y (f,tU 
dt> \ dt + dy /' 

dE' x (?,t?) _ r , dE x (f,t) | v dB z {r,t) \ 
dz' V dz c dz / 

a^(rV') ar/ _ « _ ^ ra % , . , , ^fdE z (r,t) vdB x (?,t) 



dx' dx^ z ^^ c Bx(r,t)^J ^ ^^B x (r,t) + t(^ ^ r 

^r( dV (f^ *f\_LT^ dv \( dEz vdB *\ ir( dV ut , ^\ , r^' dv\(dE z v dB x \ 

+r(-(, + ,t ) + rt - ) (— - j + r(-(t + ^) + r-- j (— - ) . (A9) 

Combining Equations IA9| and using the expressions for the x derivatives of v and T shown in Equation IA31 it is 
straightforward to obtain 

- -,cv< x ev ,n h - + »( rf f + i a -§). (mo, 

Faraday's equation: z component 

The derivation in this case is analogous to the case of the y component. From Equations (|A2|) we can get: 



B' z (f,t') = T{B z {r,t) + lE x {r,t)), E' x (f \t') = T(E x (r,t) + ^B z (r,t)), and E' y {7 ,t') = E y (f,t), 

(All) 

which imply that 

dB^.t') __ r2 / dB z (f,t) | v dE x (r,t) | v dB z (r,t) | v 2 dE x (r,t) \ 
dt' \ dt c dt dy c dy J 

dE' x (?J) _ v2( v dE x (r,t) , v 2 dB z (f,t) , dE x (r,t) , v dB z (r,t) 



/ v dE x (r,t) v 2 dB z (r,t) dE x (f,t) v dB z (r,t)\ ^ 
Vc 2 <% c 3 dt ' cty c dy J' 



dy' Vc 2 9t c 3 9t 9?/ c 9y 

«5<^> + ?ma + v() + ^ + asm ( * ( , + • „ + - 4 y {m 

Then, combining Equations (|A12|) with the derivatives of v and T given in IA31 it can be shown that 

-|3.. (( v,i V/))! . s( ^ + ^). 

The three components of Faraday's equation (Equations IA6| [A9| and lA12|) can be combined in a single expression: 

d&ifj) A A , / ,dE'(?,t') y' dE / (r f ,t')\ A /A _ N 
— = — cV x E (f ,t ) — sB x {f ,t )y + s^et + ^ x ( A 14) 

Equation IA14I contains a time dependent term arising from the time evolution of the S' coordinates with respect to 
the ones in S. Indeed, as time goes on, dE' /dx' (dB' /dx') must increase its difference relative to dE/dx (dB/dx) 
simply because two points of equal y in S (Ay = 0) have a difference in y' (Ay') that grows linearly with time. This 
explains why the time dependent term must be proportional to the extent to which the fields depend on y' . Thus, in 
the 2D case (relevant to this work), this time dependent term does not play any role. 

In addition to the time d epend ent term, there is the term proportional to y' , which appears to be inconsistent with a 
2D treatment of equation IA14I However, in the small box limit (i.e., when y' is much smaller than the distance from 
the box origin to the disk center, ro), the magnitude of the velocity sy' is much smaller than the orbital velocity of 
the disk at ro (^o)- Thus, since vo <C c, sy' must also be <C c, which allows us to neglect the space dependent term 
in Equation IA14[ especially if one expects \E'\ <C \B'\ in non-relativistic t urbu lence. Thus, given our non-relativistic, 
small box approximation, we can neglect the y' dependence in Equation IA14[ which allows us to formally use a 2D 
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approach to this problem^ In that limit, Faraday's equation can be expressed as: 



Bt' 



-cV x E'{?j) - sB' x {?j)y. (A15) 



Ampere's equation: x component 

We want to get an equation analogous to 

BE (r t) 

g ' 1 = (cV x B(f, t)) x - 4irJ x (A16) 

in the coordinate system S'. From Equations IA2I we get 

E' x (r*,t') = T(E x (r,t) + z Bz (r,t)), B' z {? \t') = T(B z (r,t) + *E x {r,t)), and B' y (j*,t') = B y (r,t). 

«-L (A17) 

Thus, from Equations (|A1|) we get that 



dE' x (fJ) _ 2 dE x (f,t) v dB z {r,t) dE x {r,t) v 2 dB z (f,t) 

dt' [ dt c dt dy c dy h 

dB'yjr 1 ,t') _dB y (r,t) 



dz' dz 



and 



dB' z {fJ) 2 v dB z (r,t) v 2 dE x {r,t) 3B z (r,t) v dE x {f,t) 

By' V Bt + c 3 ^ + + c cfy j ' 1 j 

Then, combining Equations (|A19|) and assuming that x— current in S f transforms as J x = J Xl it is possible to show 
that 

BE' ir* f) 

% ' = (cV x B>{7,t>)) x - 4irJ>. (A19) 
The assumption J' x = J x will be checked below by considering the way charge density transforms under Equations I All 

Ampere's equation: y component 

Here we proc eed similarly as for the x component of the Ampere's equation. From the field transformations defined 
in Equations (jA2|) we get that 



E' y (f,f) = E y (r,t), B' z (f,t') = T(B z (r,t) + ^E x (f,t)), and B'^f ,t') = T(B x (f,t) - ^E z (r,t)). 

(A20) 

Thus, from Equations IA1I it is possible to show that 

dE> y {r>,t>) ( dE y (f,t) dE y (f,t) } 
dt' V dt + dy r 



8B' X g , t') _ r / 6B X (r , t) v 8E Z (r, t) - 
B z^ \ Bz c Bz ' 

dB' z {fJ) ,dB z {r,t) , vdE x (r,t) 



Bx' V Bx c Bx 

fdT . . _ ,dv\ (BB Z v BE X \ {BY,. vy\ Ty'Bv\/BB z v8E x \ , , ^ 

+ (d~Jy + ^ + n + cm) + + 1> + ("ST + c ' (A21) 

Combining Equations IA22I and assuming J' y = T(J y — vp c ) (which we will check below) we obtain 

= (cV x BV,t')) y - 4nJ' y - sE<(7,t') ~ s(ct'^ + t*£). (A22) 

Ampere's equation: z component 

Finally, from the transformations lA2l we get: 

££(f,f) = r(£k(?,t) B;(t*,t') = T(B x (f,t)-^E z (r,t)), and fl^f) = B y (f,t). 

(A23) 

5 Using a Galilean instead of Lorentzian transformation of co- rr-TTl i + u • i n ,i r i-, r 

. ii2ii *j_ i . -i , , i / , , . equation IA14I does not appear, which confirms the validity of our 

ordmates and fields, it easy to snow that the y dependence m T+ . ~r — i +• • +• r •+ 

' J & r result m the non-relativistic limit. 
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Also, from I All we obtain 

dE' z (f,f) 2 dE z (f,t) v dB x (f,t) dE z (f,t) v 2 dB x (f,t) 

dt' { dt c dt dy c dy >' 

dB' x (r*, t') =p2 v dB x (r, t) v 2 dE z (f, t) | dB x (f, t) v dE z (f, t) ^ 

dy' c 2 dt c 3 dt dy c dy 

- ^ ♦ -) + o t + (fc <-) 

Combining Equations (jA25|) and assuming J' z = J z (which will be checked below) it is possible to show that 

«^) = ( ^ x * ( ,,0).-^ + .(^ + ^). 

Thus, combining the three components of Ampere's law (Equations IA191 IA221 and IA25) ) we get 

dE'(r',t') „. ~ n , , ? . ... / .dB'f?,?) y' dB'(r',t')\ A /A _ N 
= cV ' x S'^^') - 4ttJ' - sE' x {f,t')y - s(ct' + ^ x x. (A26) 

As with Faraday's equation, Ampere's equation also gets a time dependent term that disappears under the 2D 
approximation. Similarly, the 2D limit requires the y' term to be negligible so that the equations governing the 
evolution of E r are independent of y r . As in the case of Faraday's equation, this can be done by noticing that, 
since dB'/dt' w V' x E f , then V'x5'» sy'dB'/dt', provided that sy' <C c and \E'\ <C \B'\. This means that the 
y' — dependent term sy'dB'/dt' can also be neglected in the case of Ampere's equation. It is important to notice, 
however, that sy'dB'/dt' is not necessarily much smaller than dE' jdt' . Thus, although neglecting the y' — dependent 
term does not modify the MHD-scale dynamics (e.g., that of the MRI), it does not fully capture the evolution 
of the curl-less electric field component. Thus, our approach neglects the appearance of electric charges that in 
principle could influence the microphysics of the plasma. These charges, however, appear to be smaller than the ones 
due to frame rotation (see discussion in ^2j), which we are already neglecting in the context of the kinetic MRI dynamics. 

Thus, assuming a 2D geometry, Ampere's equation becomes: 

dE '% ,if) = cV x B'(7,t>) - Air? - sE' x (f,t')y. (A27) 

Particle Momentum Evolution In The S' Frame. 

The modified version of Ampere's equation assumes that currents transform as 

J' x = Jx, Jy = T(J y -vp c ), and J' z = J z . (A28) 

Based on the coordinate transformation [All it is possible to show that the charge density in S f (p' c ) is related to the 
one in S by 

p' c = T(p c - vJ y /c 2 + sy' J,/(c 2 r)), (A29) 

where u is the fluid velocity in S. Equation [A29] can be simplified under the assumptions v, sy' <C c, implying that, in the 
non-relativistic limit, charge densities do not transform. This means that simple Galil ean v elocity transformations are 
enough to reproduce the non-relativistic version of current transformation equations IA28i However, since individual 
particles in our simulations are allowed to become relativistic, we will use the relativist ic transformation of their 
momenta, p. This way, we ensure that no particles in S' exceed the speed of light. Thus, 

Px = Vx, p'y = T(p y -vy), , A3Q , 

p' z = p z , and 7 7 = r(7 - vp y /{mc 2 )), ^ ' 

where m is the mass of the particle. These relations (along with the time transformation shown in Equation IA1|) imply 
that 

dp' x dp x jdt 



dt' T(l - vuy/c 2 ) + su x y'/c 2 ' 



dp'y dpy/dt — mvdj/dt + sp x (l — vuy/c 2 ) 
W = 1 - vuy/c 2 + su x y f /(c 2 T) ' an 

dp' z dp z jdt 
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where u x and u y are the x and y velocities of th e particle in the S frame. It is straightforward to show that, in the 
non-relativistic limit (v, sy' < 1), Equations I A31I correspond to the standard force transformation between two inertial 
frames, plus an extra force sp' x y. Thus, given that the electric and magnetic fields transform analogous to the usual 
relativistic field transformation, one can easily show that the time evolution of the particles momenta in S' will be 
given by the Lorentz forces, plus the extra term sp' x y. 

In addition to these forces, we also need to find out what is the transformation of the coriolis and tidal forces forces 
on the particles. We know that in the S fram43 

du 

— 3ojqXx - 2cJ x u. (A32) 

It is straightforward to show that, in the cold limit, the S' version of Equation IA32I becomes 

du' / a 

— = -2ujqz x u. (A33) 

Since this force doe s not do work on a particle, it can be equivalently expressed in terms of pf instead of u! . Thus, 
combining Equation IA33I with the transformation to Lorentz forces found above, we get 

drf 1 it' 

= 2co oP ' y x - -coop'J + q(E' + - x B'), (A34) 

where q is the charge of the particle. 

Finally, it is important to point out that the momentum transformation shown in Equation s IA3QI is not entirely 
consistent with the evolution of particles positions. Indeed, by directly differentiating Equation [AT] one gets 





r(i 


- VUy/C 2 ) 


+ su x y' /c 2 ' 




Uy - V + 


su x t' yr 


(i- 


VUy/C 2 ) ~{ 


- su x y'/(c*T) 




u z 




r(i 


- vu v/ c2 ) 


+ su x y'/c 2 ' 



and 

(A35) 

In the non-relativistic regime, the discrepancy appears only in u' y . In the 2D case, this does not imply an inconsistency 
between the values of u' y and the evolution of the y' position of particles (where by definition y' = 0). In 3D, on 
the other hand, this implies a violation of charge conservation, which will need to be taken into account in future 3D 
studies of this problem. 



This expression is valid m the cold limit, i.e |n| << \v \ Since fluid sense Th as { as the mean velocities of the part i c l es is 
in our simulationsindividual particles are allowed to get accelerated non -relativistic, the contribution of gravity to particles dynamics 

to velocities \u\ » | Vo |, the validity of the cold limit will be m a ^ be ^ described b Equation |A32] 



